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We study the interface between a solid trapped within a bath of liquid by a suitably shaped 
non-uniform external potential. Such a potential may be constructed using lasers, external electric 
or magnetic fields or a surface template. We study a two dimensional case where a thin strip of 
solid, created in this way, is surrounded on either side by a bath of liquid with which it can easily 
exchange particles. Since height fluctuations of the interface cost energy, this interface is constrained 
to remain flat at all length scales. However, when such a solid is stressed by altering the depth of 
the potential; beyond a certain limit, it responds by relieving stress by novel interfacial fluctuations 
which involve addition or deletion of entire lattice layers of the crystal. This "layering" transition is 
a generic feature of the system regardless of the details of the interaction potential. We show how 
such interfacial fluctuations influence mass, momentum and energy transport across the interface. 
Tiny momentum impulses produce weak shock waves which travel through the interface and cause 
the spallation of crystal layers into the liquid. Kinetic and energetic constraints prevent spallation 
of partial layers from the crystal, a fact which may be of some practical use. We also study heat 
transport through the liquid-solid interface and obtain the resistances in liquid, solid and interfacial 
regions (Kapitza resistance) as the solid undergoes such layering transitions. Heat conduction, which 
shows strong signatures of the structural transformations, can be understood using a free volume 
calculation. 

PACS numbers: 68.65.-k, 68.08.-p, 62.50.+p, 44.15. +a 



I. INTRODUCTION 



Interfaces between co-existing phases in condensed 
matter systems posses an intrinsic width determined 
by competition between bulk and surface energies [1- 
4] which is usually further broadened by capillary 
fluctuations[5, 6] viz. random fluctuations of the inter- 
face away from its mean position. The latter is a direct 
manifestation of the fact that translations of the interface 
in a direction perpendicular to its plane cost no energy [7]. 
Such capillary fluctuations actually cause the interfacial 
width to diverge [7, 8] with a cutoff determined by the sys- 
tem size. External potentials eg. gravitational helds or 
coupling to a substrate however [9-11] tend to round-off 
this divergence and set additional length scales, absent 
in the free system. If one or both of the phases hap- 
pens to be a solid, then long-ranged elastic energy costs 
for deforming the interface may also limit the broad- 
ening of the interface in equilibrium [12]. Nevertheless, 
the analog of capillary fluctuations, viz. crystallization 
waves have been experimentally observed at liquid-solid 
interfaces [13]. It is, of course, also possible to suppress 
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interfacial capillary fluctuations by position dependent 
chemical potentials which break the translational sym- 
metry of the interface explicitly [14]. One would expect 
that such interfaces would remain essentially flat over all 
length scales and therefore be completely inert. 

We show, in this paper, that this is not so. Surpris- 
ingly, there are novel fluctuations and phenomena asso- 
ciated with such constrained interfaces which have static 
as well as dynamic consequences. Particles are trans- 
ferred across the interface in new and interesting ways. 
A liquid-solid interface constructed in such a fashion can 
have physically interesting height fluctuations while at 
the same time remaining flat! The system manages this 
by selecting only those fluctuations which involve changes 
in height by a single atomic spacing and have a wave- 
length equal to the size of system in the transverse direc- 
tion. Further, such interfacial fluctuations are driven by 
the elastic response of the solid to stresses imposed by 
the external potential. As the depth of the trapping po- 
tential is gradually increased, the solid accomodates this 
stress by incorporating layers of atoms from the liquid, 
either only parallel to the interface or alternately perpen- 
dicular and parallel in a cyclical fashion, depending on 
the nature of the interactions. We believe that some of 
our predictions may be directly checked for liquid-solid 
interfaces in atomic, as well as, colloidal systems where 
the chemical potential field may be provided either by a 
laser trap or by a patterned substrate. Needless to say, 



such fluctuations are expected to be observed only in sys- 
tems where the overall size is small - of the order of only 
a few atomic spacings in the transverse direction. Pre- 
liminary results of our studies of this system have been 
published in Refs.[15] and [16]. 

This paper is organized as follows. We begin, in 
the next section, by showing how a two-dimensional 
liquid-solid interface may be produced using a non- 
uniform external potential. Monte Carlo (MC) com- 
puter simulations [17] in the constant number, area and 
temperature (NAT) ensemble are set up to realize this 
explicitly for particles interacting with hard disk, soft 
disk and Lennard- Jones potentials. The interface is 
then characterized using a variety of thermodynamic 
and structural quantities, which are measured as a func- 
tion of the perpendicular distance from the interface. 
As a function of the depth of the potential well, the 
trapped solid undergoes, what we have called, "layering" 
transitions[15, 16, 18-21] which involve the addition (or 
removal) of an entire layer of solid from (or into) the sur- 
rounding liquid through the interface. This transition is 
described in detail for the hard disk system in section III. 
The layering transition is accompanied by a sharp jump 
in the density of the solid. We obtain this density jump 
within a mean field, thermodynamic approach[20]. A 
comparison of the predictions of the thermodynamic the- 
ory and our MC computer experiments show the theory 
to be approximately correct to within a few percent. We 
show that the layering transition is a novel mechanism by 
which a stressed nano-solid constrained by an external 
potential can respond plastically to large stresses with- 
out nucleating dislocations or cracks [22, 23] . We establish 
that this phenomenon is general and is independent of the 
particular interatomic potential used. In section IV, we 
describe the kinetics of the layering transition. In section 
V we explore how mass and momentum are transported 
across the liquid-solid interface and especially the role 
of the layering transition on the transport coefficients. 
Molecular dynamics (MD) simulations in the constant 
number, area and energy (NAE) ensemble[17] are car- 
ried out for this purpose. We show in this paper that 
(1) fluctuations associated with these transitions are of 
a special kind always involving the transfer of complete 
layers of solid (2) these fluctuations offer resistance to 
the transfer of momentum and energy through the in- 
terface (3) the resistance is maximum when the energy 
matches that required to raise a complete lattice layer 
from within the potential well into the surrounding liq- 
uid. We study the stability of surface kinks at the liquid 
solid interface in the hard disk system. We then study 
the response of the interface to weak acoustic shocks [24] 
which are shown to cause spallation of complete lattice 
layers if the incident energy of the shock is large enough. 
In section VI we use non-equilibrium molecular dynam- 
ics to study heat transport [21] through the liquid-solid 
interface in soft disks and obtain the conductance of the 
system in different regions, heat current and contact or 
Kapitza resistance[25, 26] of the interface as a function 
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FIG. 1: (Color online) A schematic diagram of the system 
showing the liquid and solid regions produced by the exter- 
nal chemical potential of depth — /j,. The various dimensions 
mentioned in the text are also marked in the figure. 



of the depth of the potential well. The heat conduction 
is particularly sensitive to the fluctuations in the direc- 
tion of current flow. We present an approximate free 
volume type calculation that qualitatively captures the 
response of the heat conductance to the internal changes 
in the solid induced by the trapping potential[21]. Fi- 
nally, in section VII we conclude after discussing some 
consequences of our study and its relevance to experi- 
ments. 



II. CONSTRAINED LIQUID-SOLID 
INTERFACE: STATIC PROPERTIES 

In this section, we explore the possibility of creating 
a patterned sequence of confined solid and liquid re- 
gions using an external, space-dependent, chemical po- 
tential field </>(r). Consider a two dimensional system (see 
Fig.l) of N atoms of average density (packing fraction) 
i] = ttN/4A within a rectangular cell of size A = L x x L y 
where the central region S of area A s — L x x L s is oc- 
cupied by N s atoms arranged as a crystalline solid of 
density r/ s > 77, while the rest of the cell is filled with 
liquid of density r\g < r\. The difference in density is 
produced by an external field 0(r) = — /j, for r £ S; in- 
creasing sharply but smoothly to zero elsewhere with a 
hyperbolic tangent profile of width 8^. How may (j)(r) 
be realized in practice? In model solids like colloids [27], 
one may use a surface template to create a static pattern 
[28]. In real systems, as well as colloids [29], one may 
be able to use laser traps [30] or non -uniform electric or 
magnetic fields. Usual laser traps for alkali metals or rare 
gas atoms are in the range of lOmK for which a power 
of about lOOmW is required[31]. For colloidal systems, 
required laser powers are even lower [3 2]. 

We first describe our results for a system where the 
atoms interact with a hard disk potential [33] which is 
infinitely repulsive if the distance r.y between two atoms 
i and j is less than a, the hard disk diameter and zero 
otherwise. We show later that qualitative results for more 
realistic potentials [17], e.g., soft disk or Lennard- Jones 
are similar. We have chosen 8$ = o-/A, where a is the 
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FIG. 2: The density profile rj(y) coarse grained over strips 
of width a (averages taken over 10 3 MC configurations each 
separated by 10 3 MCS) at p — 6. The dotted lines show the 
predictions for liquid density r\i — 0.6725 and solid density 
rj s — 0.752 from a simple thermodynamic theory presented in 
Sec.III.A. 




FIG. 3: (Color online) Solid-liquid interface at p = 8. Super- 
position of 500 configurations separated by 10 3 MCS showing 
a solid like order (red : high rj) gradually vanishing into the 
fluid (blue : low rj) across a well defined solid-fluid interface. 



hard disk diameter and sets the scale of length. The 
energy scale for this system is set by ksT where k B is 
the Boltzmann constant and T the temperature. In our 
simulations we set a — 1 and fc^T = 1 unless otherwise 
stated. 

The full configuration dependent Hamiltonian is Ti — 
Ylij u ( r ij) + Si^( r i)- We have carried out extensive 
MC simulations with usual Metropolis moves [17], peri- 
odic boundary conditions in both directions and in the 
constant number, area and temperature ensemble to ob- 
tain the equilibrium behaviour of this system for differ- 
ent /i at a fixed average r). N = 1200 particles occupy 
an area A = 22.78 x 59.18 with the solid occupying the 
central third of the cell of size L s = 19.73. The initial 
configuration is chosen to be a liquid with rj = 0.699; 
close to but slightly lower than the bulk freezing density 
r]f = 0.706 [33]. On equilibration, S contains a solid 
with the close-packed planes parallel to the solid-liquid 
interfaces which he, at all times, along the lines where 
4>(y) — > 0. The equilibration time is large and many 
(~ 10 7 ) Monte Carlo steps (MCS) are discarded before 
results shown in Figs. 2 to 7 are obtained. 

The density rj(y), coarse grained over strips of width 
~ cr, varies from its value rji in the liquid to rj s as we 
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FIG. 4: Bond orientational order parameter across the liquid- 
solid interface for a 21 layered solid at p — 6 surrounded by 
liquid on both sides. 

move into the region S (Fig. 2). Averages are taken over 
10 3 MC configurations each separated by 10 3 MC steps. 
The trap depth p = 6, supports an equilibrium solid 
of density ?y s = 0.753 in contact with a fluid of den- 
sity Tjg = 0.672. The horizontal lines are predictions of 
a simple free- volume based theory ([15]) for r/ s and rj£. 
We discuss the theory in Sec. III. A. A superposition of 
atomic positions shows a static, flat, liquid-solid inter- 
face with the solid like order gradually vanishing into the 
liquid (Fig. 3). We have thus created a thin nano-sized 
crystal which is 21 atomic layers wide (for a trap depth, 
p, = 6) and is flanked on either side by liquid separated 
by two liquid-solid interfaces. 

The bond orientational order parameter (ipe(y)) where, 
the local value of ipe for a particle i located at = (x, y) 
is given by 

1 3 

The sum is over the j € Ni neighbours of i-th particle, 
and 8ij is the angle between the vector ry and an ar- 
bitrary but fixed reference axis. To obtain (ipgly)), this 
quantity is coarse grained over strips of width a (aver- 
ages taken over 10 3 MC configurations each separated by 
10 3 MCS). This shows a sharp rise from zero to a value 
close to one, as we move into the region S (Fig. 4). This 
indicates that the particles in S maintains hexatic order. 
However, this does not necessarily justify the phase to be 
a solid. Therefore, there is the need to calculate the solid 
order parameter. 

The order parameters corresponding to the solid phase 
are the Fourier components of the (nonuniform) density- 
density correlation (p(ri)p(rj)) calculated at the recip- 
rocal lattice points {G}. This (infinite) set of numbers 
are all zero (for G =/= 0) in a uniform liquid phase and 
nonzero in a solid. We restrict ourselves to the star con- 
sisting of the six smallest reciprocal lattice vectors of the 
two dimensional triangular lattice. In modulated liquid 
phase, the Fourier components corresponding to two out 
of these six vectors, eg., those in the direction perpen- 
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FIG. 5: (a)The reciprocal lattice vectors Gi and G2, and the 
rectangular unit cell, (b) Solid order parameter corresponding 
to Gi is nonzero in the region S at p = 6, indicating a solid. 
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FIG. 6: From left to right, the three view-graphs show struc- 
ture factors for the liquid, the solid and the interfacial regions 
respectively at p = 6. 



dicular to the interface, G2, arc nonzero. The other four 
components of this set which are equivalent by symmetry 
(Gi) are zero in the (modulated) liquid and nonzero in 
the solid (if there is true long ranged order). Thus, we 
use the following order parameter : 
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I cxpHG fc • ry)|) 



(2) 



where = Vi — Vj. The solid order parameter so defined, 
in the direction Gi (see Fig. 5) and the others equivalent 
to it by symmetry is nonzero in the region S, indicating 
the nucleation of a solid phase. Note that {xpo) shows a 
larger interfacial region than that obtained from (f^Gi)- 
This is because the liquid near the liquid solid interface 
is orientationally highly ordered due to the proximity of 
the solid. 

The structure factor describes density correlations in 
Fourier space, 



S(k) = N- 1 (p(k)p(-k)) 



(3) 



where p(k) = ex P(*k ' r «)- ^ n a simulation with 

periodic boundaries, k is restricted by the periodicity 
of the system, i.e. with the simulation box. The two- 
dimensional structure factor shows sharp peaks at trian- 
gular lattice positions for the solid region and isotropic 
ring pattern for the liquid (Fig. 6). The interfacial region 
also shows diffuse peaks at approximately triangular lat- 
tice positions indicating once more that the interfacial 
region has considerable orientational and short ranged 
translational order. 

Before we end this section, we must voice a note of cau- 
tion about the identity of phases in small and confined 



systems similar to the one we have here. It is well known, 
for example, that in two dimensions it is impossible to 
obtain a solid with true long ranged order [34]. Displace- 
ment correlations defined as (u(r) • u(0)), where the dis- 
placement vector u is measured with respect to the zero 
temperature perfectly crystalline reference solid, grow as 
log(r). This implies that true Bragg peaks are impossi- 
ble in two dimensions. Nevertheless, finite size and lack 
of complete averaging can lead to structure factors with 
sharp peaks, while obtaining the true logarithmic diver- 
gence may need considerable amount of computational 
effort. For a solid confined to a two-dimensional channel, 
the situation is even more dramatic. The displacement 
correlations now increase linearly with system size [35] for 
distances larger than a crossover length ~ L s log(L s ). 
The structure factor should show true Bragg peaks for 
reflections from planes parallel to the confining walls and 
should be diffuse in the other directions, implying there- 
fore smectic like ordering throughout. In this paper, how- 
ever, we continue to use the words "solid" and "liquid" in 
the usual sense referring to the presence or lack of solid- 
like order as shown in Fig. 3 and Fig. 5. This is mainly 
due to the fact that our motivation of this study is to 
probe and understand the properties of small (nano) sys- 
tems. For such systems with channel length L x not much 
larger than L s the phonon fluctuations are suppressed 
and can not destabilize the solid. However, even for very 
large strips where phonon fluctuation destabilizes a solid 
to a smectic phase, layering transitions are expected to 
occur [36] and some associated proterties like accoustic 
spallation of layers and large change in heat conduction, 
which we describe in this paper, are expected to remain 
operative. We must, at the same time, keep in mind that 
some of the properties of the confined solid, including the 
layering transition and the case of spallation of solid lay- 
ers in response to weak acoustic shock (to be discussed 
below), may in- fact, be a consequence of incomplete or- 
dering. 



III. THE LAYERING TRANSITION 

We now calculate the difference in densities between 
the solid and liquid regions A 77 = (i] s — i]e) as a function of 
the strength of the external field p. While A77 jr\ increases 
with increasing p as expected, the smooth increase is 
punctuated by a sharp jump (Fig. 7). An examination 
of the particle configuration shows that the jump occurs 
when an extra close -packed layer enters S increasing 
the number of solid layers by one. For the parameters 
in our simulation, the jump occurs at p w 8 with the 
number of layers increasing from 21 to 22. The value 
of p at the jump is a strong function of L s . The solid 
structure is seen to be a defect free triangular lattice 
with a small rectangular distortion ed(f] s ,L s ) [20]. Wc 
have examined the variation of Ar](p)/r] by cycling p 
adiabatically around the region of the jump. This yields 
a prominent hysteresis loop as shown in Fig. 8 which 
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FIG. 7: Plot of the equilibrium fractional density change 
An/n as a function of n (points (MC data), solid line (free 
volume theory), showing discontinuous jump at fi « 8. The 
labels A-D mark the stable 21 layer solid (A), the transition 
(B-C) and the stable 22 layer solid (D) respectively. 
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FIG. 8: Hysteresis loop as p is cycled at the rate of 0.2 per 
10 6 MCS. The central jagged line is the result of the initial 
cycle when a single dislocation pair was present in the solid 
region. 



indicates that 'surface' steps (dislocation pairs) nucleated 
in the course of adding (or subtracting) a solid layer, 
have a vanishingly short lifetime. Consistent with this we 
find that the jump in An(fj,)/r] vanishes when the system 
is minimised at each fi with a constraint that the solid 
contains a single dislocation pair (Fig. 8). Interestingly, a 
dislocation pair forced initially into the bulk, rises to the 
solid- fluid interface due to a gain in strain energy [37], 
where they form surface indentations flanked by kink- 
antikink pairs. This costs energy due to the confining 
potential, as a result, the kink-antikink pair gets quickly 
annihilated by incorporating particles from the adjacent 
fluid. The jump in Arj(^)/n is also seen to decrease with 
increasing 5$. 



A. The thermodynamic theory 

The qualitative features of these results may be ob- 
tained by a simple thermodynamic theory (Fig. 7) with 



FIG. 9: In our free- volume calculation we assume that the 
outer six disks are fixed at their average positions and the 
central disk moves within this cage of fixed particles. The 
curve in bold line shows the boundary B of the free volume 
available to the central test particle. A point on this boundary 
is denoted by Po(x,y) while the centers of the six fixed disks 
are denoted by Pi(xi,yi) with i — 1,2. ..6. b and h denotes 
the base and the height which uniquely decides the perimeter 
of B and the enclosed free volume. These are functions of 
density ij s and width of the potential well L s . 



harmonic distortions of the solid, ignoring contributions 
from spatial variations of the density. Note that this the- 
ory is in the spirit of a mean field approximation where 
all effects of fluctuations discussed in the concluding part 
of the last section are ignored. The free energy of the to- 
tal system is written down as a sum over the free energies 
of the solid and the fluid. The free energies of the bulk 
hard disk fluid and solid are relatively easy to obtain as 
shown below. 



1. Free energy of the solid 

In order that the solid channel accommodates n/ layers 
of a homogeneous triangular lattice with lattice parame- 
ter ao of hard disks of diameter a, the channel width 
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Defining 



x(Vs,L s ) = 1 + 



2(L a -a) 
\/3a 



(4) 



(5) 



if x = integer — n\ Eq.4 is recovered and the channel can 
accomodated ni layers of homogeneous triangular lattice 
and x 7^ integer implies a rectangular strain away from 
the perfect trianglular lattice. For any given set of values 
for L s and n s one can find \ in the following manner. One 
can associate a triangular lattice of lattice parameter a 
from any given n s = 7r/2\/3ao- This set of ao and L s 
defines a specific X- Then the channel can accommodate 
n ; = int(x) (nearest integer to x) number of layers of a 
centered rectangular (CR) lattice with lattice parameters 
a y = 2(L S — u)/(ni — 1) and, a x — ir/2n s a y . This lattice 
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has strains e xx — ^ 
deviatoric strain Sd = £ a 



1 and e 



yy ~ 
E yy is then, 
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ni 



1 



*=lr - 1. The 

ni-l 



(6) 



In the fixed neighbour free volume theory (FNFVT), 
particles of high density solid are assumed to be con- 
fined within a cage formed by the average positions of 
its nearest neighbours. This cage and the available free 
volume of the test particle to move around in this cage 
is entirely defined by the quantities b = a (l + e xx ) and 
h = (V3a /2)(1 + e yy ) (Fig.9)[21]. The amount of avail- 
able free volume Vf v (r} s ,L s ) bounded by the bold line B 
in Fig. 9, can be calculated by using tedious but straight- 
forward geometrical considerations. The free volume free 
energy density is f s {rj s ,L s ) = -(4%/tt) k B T\og(v fv ). f s 
always remains an upper bound to the exact free energy. 
This upper bound becomes assymptotically exact in the 
close-packed limit. Since e xx and e yy have discontinuity 
at half integral x, the free energy f s has maxima at those 
X values[20]. In Sec. VI we present an approximate the- 
ory for calculating heat conductance in solid using this 
free volume approach [21]. 



2. Free energy of the liquid 

The free energy density of the liquid bulk phase may 
be simply written as 



ft = Pi 
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(7) 



where pi = ir/i/ir, the ideal gas Helmholtz free energy 
density f id = p e \og(p e ) — p e and we use the semiemperical 
equation of state for the hard disk liquid in Ref . [38] . This 
equation of state has been observed to show excellent 
agreement with hard disk fluid up to p t = 0.65 even 
when the fluid is confined in a hard narrow channel [39]. 
Following Ref. [38] we use f t = f san + f id where 



fsan Pi 



(2 Vc - 1) log (l - (2ri c - - log(l - fj 



2(1 -Vc) 



(8) 



with rj c — 7r/2\/3, the close-packed density for 2D hard 
disk triangular lattice. 



3. Free energy of the system 

We now write down the total free energy density of 
the system (fluid + solid regions) using the free energy 
density expressions for solid and liquid bulk phases as 

/ = x [f s (Vs,L s ) - 4j? g /i/7r] + (1 - x)f e (r)t)- (9) 
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FIG. 10: A plot of the deviatoric stress 7 d against strain 
Sd- The arrows show the behaviour of these quantities as 
p is increased from the points marked A - D. The labels 
correspond to the same states as in Fig. 7. 



We then minimise this free energy density with the con- 
straint that the average density is fixed, rj = xr\ s + (1 — 
x)rje, where x is the area fraction occupied by S. The 
result of this calculation is shown in Fig. 7 where it is 
seen to reproduce the jump in Arj(p)/rj. 

Why does the solid incorporate layers of atoms from 
the liquid ? This question may be answered elegantly if 
one calculates the deviatoric stress "fd in the solid region 
as a function of the depth of the strain, Ed- The stress 
may in fact be obtained in a straight forward fashion from 
the expression of the free energy. Differentiating the free 
energy of the solid with respect to s d we obtain 



Id 



ds d 



(10) 



When 7d is plotted versus the deviatoric strain e dl we 
observe that the solid is not stress free for any arbitrary 
combination of p and L s . In fact, for our parameters, 
initially the 21 layered solid is under tension in y- di- 
rection. We follow the variation of the deviatoric stress 
with the strain as p is increased from the points A — D 
in the Fig. 10. The state of stress in the solid jumps 
discontinuously from tensile to compressive from B — > C 
due to an increase in the number of solid layers by one 
accomplished by incorporating particles from the fluid. 
This transition is reversible and the system relaxes from 
a state of compression to tension by ejecting this layer as 
p is decreased. As p is increased, the tension increases till 
7d reaches about —2.96 when the corresponding strain is 
about —0.052. At this point a layer enters the solid region 
and the stress and strain switches from tensile to com- 
pressive. Further increase in p now decreases the stress 
and drives the solid to a state of zero stress at p = 10. 
Thus, the layering transition from 21 to 22 layers as ob- 
served by us is a mechanism for relieving stress. 

In our theory we assumed the channel can accomo- 
date only an integer number of layers and we ignored 
any possibility of defects. However, our theory allowed 
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FIG. 11: Plot of the equilibrium fractional density change 
Ap/p as a function of p (points - MC data for soft core; line 
is a guide to eye) . 
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FIG. 12: Plot of the equilibrium fractional density change 
Ap/p as a function of p (points - MC data for Lennard- Jones 
system; line is a guide to eye). 



for continuous change of base length b. In small solid, 
this change is far from continuous. In our FNFVT, while 
h remains constant until a layering transition, b contin- 
uously reduces with increase in p. In reality, b can only 
be reduced by incorporating a new particle in an exist- 
ing layer. Once a solid is formed (for pi > 4, see Fig. 7), 
formation of point or line defects are energetically very 
costly, thus even on an average b can reduce only if all 
the layers include a particle each and lattice parameters 
shrink in coherence. This is the reason why in theoreti- 
cal prediction Arj/rj grows smoothly in between the layer- 
ing transitions, though the same quantity remains almost 
constant in simulation (Fig. 7). We shall see in later sec- 
tions that the coherent absorption of particles by all the 
layers at a time is indeed observed at larger well depth pi 
in similar simulations with soft core particles. 



ulation e = a = 1 and N — 1200 particles occupy an 
area A = 24 x 60 with the solid occupying the central 
third of the cell of size L s = 20. The average density 
of the system is therefore kept at p — 0.833 which is to 
be compared with the freezing densities p ~ 0.987 [40] 
and p ~ 0.865 [41] for the soft core and Lennard- Jones 
systems respectively. For the soft core potential, the pi 
value at the jump in density (Fig. 11) is even quantita- 
tively comparable to the corresponding hard core system. 
In soft core systems the other possible mode of stress re- 
laxation, namely, coherent inclusion of one particle each 
in all the existing layers is observed at a larger well depth 
of p, rts 16. We shall discuss about this farther in Sec. VI 
while discussing heat transport in this soft core system. 



IV. KINETICS OF LAYERING 



B. Layering in other potentials 

If the layering transition observed in the hard disk sys- 
tem is actually a new mechanism for relieving stress in a 
thin crystal, it should be independent of the details of the 
potential. Our main results trivially extend to particles 
interacting with any form of repulsive potential, or even 
when the interactions are augmented by a short range 
attraction, provided we choose p deeper than the depth 
of the attractive potential. 

In this section we show explicitly that the layering 
transition is present in the soft core and Lcnnard-Jones 
systems. We choose ksT = 1. Once again we perform 
MC simulations in the constant NAT ensemble with peri- 
odic boundary conditions and with the external chemical 
potential pi. The relevant parameters corresponding to 
these potentials, namely e and a in u ss (r) = e (^) = 

, set the energy 
Here r = r 51 is 



ar- 12 and u LJ (r) = 4e [(^) 12 - (^) 6 
scale and the length scale respectively. 



the distance between the pair of atoms In our sim- 



The large hysteresis loops associated with the layering 
transition obtained in the last section makes it clear that 
the kinetics of this transition is slow. To study the life- 
time of the kink-antikink pairs (surface step), we resort to 
an MD simulation, using a velocity Verlet algorithm [17], 
with the unit of time given by r — yjma 2 /ksT ', where 
m(= 1) is the mass of the hard disks. Using values of m 
and a typical for atomic systems like Ar or Rb, r w lps. 
A time step of At = 7 x 10~ 5 r conserves the total energy 
to within 1 in 10 3 (at worst). 

Starting with an equilibrium configuration for the hard 
disk system, taken from our MC runs as discussed in Sec. 
II at pi — 9.6 corresponding to a 22-layer solid, we create 
a unit surface step of length I by displacing a few in- 
terfacial atoms from the solid region into the liquid and 
'quench' across the transition to p = 4.8, where a 21- 
layer solid is stable. The rest of the parameters are kept 
identical to those given in section II. We observe that 
the fluctuation thus created rapidly relaxes back and the 
surface step vanishes as the atoms are pulled back into 
the solid. We illustrate this by plotting the number of 
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FIG. 13: Plot ol the number of particles in the solid region 
N a as a function of time t (in units of r) clearly shows that 
the displaced particles are pushed back into the region S. 



hard disks within the solid region as a function of the 
MD time steps (Fig. 13). As soon as a step is created, 
the line of atoms in the portion of the solid thus exposed 
bend to fill up the gap created between the atomic layers 
and the edge of the potential trap. This generates con- 
siderable local elastic stress. Also, the liquid layer lying 
immediately adjacent to the solid has a lot of orienta- 
tional and solid-like order. For short times it responds 
elastically to the presence of an increased local density of 
atoms. The combined elastic response therefore pushes 
the displaced atoms back into the solid region thereby 
annihilating the step. This annihilation is a transient 
response and happens within a time of t, whereas the 
lifetime of the metastable 22-layered solid is about 10r, 
which is one order of magnitude larger! Indeed, a free en- 
ergy audit involving a bulk free energy gain AF ~ 1/L S , 
going from a 22 to 21 layered solid, and an elastic energy 
cost ~ log(i) for creating a step of size I, reveals that a 
surface step is stable only if I > I* ~ 1/L S . For small 
L s , the critical size I* may therefore exceed L x , the total 
length of the interface. Of course, if the step spans the 
entire length of the interface, there is no bending of the 
atomic lines and there is no elastic energy cost. This ex- 
plains the slow kinetics since the system has to wait till a 
rare random fluctuation, which displaces all the atoms in 
a solid layer across the interface coherently, gives rise to 
layering transition. Although we have explicitly demon- 
strated this for the hard disk system, we believe that 
similar considerations should be appropriate for the soft 
disk and Lennard Jones systems too. 

The slow kinetics of the layering transition may have 
an impact on the transfer of momentum across the liquid- 
solid interface in the form of regular sound waves or 
acoustic shocks. The large effective compressibility of the 
solid at the layering transition as evidenced by the jump 
in the density as the chemical potential is increased by 
an infinitesimal amount (Fig. 8) should reduce the veloc- 
ity of sound considerably. The propagation and scatter- 
ing of sound in an inhomogeneous region with coexisting 



phases has been studied extensively [42-45] in the past. 
The transfer of mass between coexisting phases at inter- 
phase boundaries is known to slow down and dissipate 
sound waves traveling through the system. Our system 
has an artificially created inhomogeneity, which should 
have a similar effect on its acoustic properties. 

Further, the mechanism of stress relaxation of a thin 
[L a small) solid via the transfer of an entire layer of atoms 
may be exploited in a variety of practical applications, 
provided we can eject this layer of atoms deep into the ad- 
joining fluid and enhance its lifetime. We may be able to 
use the ejected layer of atoms to create monolayer atomic 
films or coatings [46] . Highly stressed mono-atomic layers 
tend to disintegrate or curl up [47] as they separate off 
from the parent crystal. It may be possible to bypass 
this eventuality, if the time scale of separation is made 
much smaller than the lifetime of the layer. Can acous- 
tic spallation [24] be used to cleave atomic layers from a 
metastable, stressed nanocrystal? 

In the next section, therefore, we study the response 
of the liquid solid interface in our system of hard disks 
to acoustic shocks with a view to studying the effect of 
the layering transition on acoustic shock propagation and 
dissipation as well as the properties of the ejected layer. 



V. MASS AND MOMENTUM TRANSFER 

Consider sending in a sharp laser (or ultrasonic) pulse, 
producing a momentum impulse (v y (t = 0) = Vo) over 
a thin region in y spanning the length L x of the sim- 
ulation cell. This results in a weak acoustic shock [24] 
(corresponding to a laser power « 10 2 mW and a pulse 
duration lps for a typical atomic system). The initial 
momentum pulse travels through the solid and emerges 
at the far end (Fig. 14) as a broadened Gaussian. The 
width of this Gaussian pulse, A, is a measure of absorp- 
tion of the acoustic energy of the pulse due to combined 
dissipation in the liquid, the solid and at the interfaces 
[42, 43, 48]. For large enough pulse strengths Vo, this is 
accompanied by a coherent ejection of the (single) outer 
layer of atoms into the fluid. Note that in our MD sim- 
ulations, to reduce interference from the reflected pulse 
through periodic boundary conditions, we increase the 
fluid regions on either side, so that for the MD calcu- 
lations we have a cell of size 22.78 x 186.98 comprising 
3600 particles. This is accomplished by separately equi- 
librating two liquid regions of appropriate size and den- 
sity and smoothly sandwiching our equilibrated system 
(which includes the solid strip) in between these liquid 
regions. The whole 3600 particle system is equilibrated 
for a further 10 4 MCS before it is used as an initial con- 
dition for the MD simulations. In Fig. 14 we show the 
initial momentum pulse with strength Vo — 6. as pro- 
duced within a narrow strip of size ~ a, just to the left 
of the solid region and the curves are fitted to a Gaussian 
(and the width A 2 extracted) when the maximum of the 
pulse reaches a fixed distance of 44.1 from the source. A 
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FIG. 14: (a)-(c) Plot of the absolute value of the momentum 
|t>j/(j/)| for molecular dynamics times t — 0.0007 (a), 0.2828 
(b) and 2.8284 (c). The dotted lines show the position of the 
solid-fluid interfaces. The fit to a Gaussian (thick solid line) 
is also shown in (c). Curves such as in (a)-(c) are obtained 
by averaging over 100 — 300 separate runs using different re- 
alizations of the initial momentum. 




FIG. 15: Plot of the squared width A 2 of the momentum 
pulse after it emerges from the solid as a function of Vo for 
fi = 4.8 ( O) and 9.6 ( □). The solid lines are fits to an effective 
liquid theory. The peak in A 2 (Vb) so produced is more promi- 
nent for the metastable 22 layered solid /i = 4.8 than for the 
stable (fi = 9.6) system showing a more coherent momentum 
transfer in the former case. 



reflected pulse can also be seen. 

When a shock wave, which propagates through a con- 
ventional solid emerges from the free surface, the com- 
pressed material expands - or unloads - to zero pressure 
[24] . The unloading (rarefaction) wave travels backwards 
into the material with the speed of sound. The response 
of the solid depends on the specific nature of the shock 
front. For a shock wave with an approximately Gaussian 
profile as in our case, significant negative pressures can 
develop at the interface where the shock emerges due to 
the interaction of the forward and the reflected waves and 
a portion of the solid may split off by a process known 
as "spallation". Spallation in bulk solids like steel needs 
acoustic pressures in excess of 10 5 N/cm 2 [24] usually 
available only during impulsive loading conditions; the 
ejected layer is a "chunk" of the surface. In contrast, 
the pressures generated by the shock wave in our system 
causing coherence nano spallation involves much smaller 
surface stresses of the order ksT/a 2 s=s 10~ 5 N/cm 2 . This 
difference comes about because unlike a bulk system, a 
strained nanocrystal on the verge of a transition from 
a metastable n + 1 to a n layered state readily absorbs 
kinetic energy from the pulse. Also, as mentioned be- 
fore, a confined solid strip in two-dimensions has very 
strong smectic ordering which effectively decreases the 
coupling between the layers[20, 35]. In fact, a quasi one- 
dimensional solid (L x S> L s ) is better regarded as a smec- 
tic with weak solid like modulations. The fact that sur- 
face indentations are unstable (Fig. 16) unless of a size 



comparable to the length of the crystal, L x , ensures that 
a full atomic layer is evicted almost always, leading to co- 
herent absorption of the pulse energy. The coherence of 
this absorption mechanism is markedly evident in a plot 
of A 2 against Vq which shows a sharp peak (Fig. 15). 
Among the two systems studied by us, viz., a metastable 
(/i = 4.8) and a stable (n = 9.6) 22 layered solid, the for- 
mer shows a sharper resonance. Note that the absorption 
of momentum is largest when the available kinetic energy 
of the pulse exactly matches the potential energy required 
to eject a layer. To elucidate this fact further, we plot the 
configurations of the metastable system (p, = 4.8) as the 
pulse travels through the system, for two different pulse 
strengths, Vq = 2 and Vo = 6 (Fig. 16). The weaker 
momentum pulse (Vo = 2) initially ejects a few atoms 
of the interfacial crystalline layer of the metastable 22 
layered solid. However, the resulting large non-uniform 
clastic strain evidenced by the bending of lattice layers 
causes these atoms to be subsequently pulled back into 
the solid. This effect is the same as that seen in last sec- 
tion. Only a stronger pulse, Vo = 6, capable of ejecting 
a complete lattice layer succeeds in reducing the number 
of solid layers by one leading to an overall lower elastic 
energy. 

The eviction of the atomic layer is therefore assisted by 
the strain induced interlayer transition and metastability 
of the 22 layered solid discussed above. Spallation is also 
facilitated if the atomic interactions are anisotropic so 
that attraction within layers is stronger than between lay- 
ers (eg. graphite and layered oxides [47]), for our model 
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FIG. 16: (Color online) (a) Configuration snapshot from a 
portion of our MD cell showing hard disk atoms (green circles) 
at the solid (bottom)- liquid interface (yellow line) as a weak 
momentum pulse (Vo = 2) emerges into the liquid, at three 
different times, for /i = 4.8. The pulse initially ejects a few 
atoms from the solid but are subsequently pulled back due 
to the large elastic strain cost in bending of the interfacial 
crystalline layer (red circles), (b) The same for a stronger 
momentum pulse, Vo = 6. This time the pulse strength is 
sufficient to eject the layer. 

of purely repulsive hard disk solid, an effective intralayer 
attractive potential of mean force is induced by the ex- 
ternal potential[20]. 

The spallated solid layer emerges from the solid surface 
into the fluid, and travels a distance close to the mean 
free path; whereupon it disintegrates due to viscous dis- 
sipation (Fig. 17). A simple estimation of the life-time 
of the ejected layer may be undertaken as follows. To 
obtain the life time of the spallated layer we obtain the 
time development of the Fourier component of the local 
density correlation (ipG 2 {y)) which is just a time depen- 
dent generalization of the quantity shown in Fig. 5. We 
obtain this by averaging, at each time slice t, the quan- 
tity Y^u j) exp(— iG 2 • (rj — rj)) over all nearest neigh- 
bour pairs with centre of the vector (r; — rj) lieing 
within a strip of width ~ a centered about y and span- 
ning the system in x. The wavenumber G2 = (2ir/d)n 
where d = .92 is the distance between crystal lines in 
the direction n normal to the fluid-solid interface. The 
solid (central region with ipG k (y,t) 7^ 0, for k — 1,2,3) 
ejects a layer (shown by an arrow in Fig. 17) which subse- 
quently dissolves in the fluid. The lifetime of the layer is 
around 2-3 time units (r) which corresponds to a few ps 
for typical atomic systems. The lifetime increases with 
decreasing viscosity of the surrounding fluid. Using the 
Enskog approximation [49] to the hard disk viscosity, one 
can calculate the bulk viscosity for a hard disk fluid to 
be [49, 50] 

(e = -Coo?? 2 .g(a) (11) 
where Coo is a constant and g(a) is the pair-correlation 
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FIG. 17: (Color online) A plot of the time development 
{ipG 2 (y))- The solid ejects a layer (shown by an arrow) which 
subsequently dissolves in the fluid. The curves from bottom to 
top correspond to time slices at intervals of At — .07r starting 
from t = 1.06t (bottom). We have shifted each curve upward 
by .03i/At for better visibility. 



function at contact. For a system of hard disks with 
m = a = (3 = 1 Coo = l/20r[5O], Thus, C,e oc rj 2 
and we estimate that by lowering the fluid density one 
may increase the lifetime of the layer considerably. The 
lifetime enhancement is even greater if the fluid in contact 
is a low density gas (when the interparticle potential has 
an attractive part [51]). 



A. Effective liquid theory 

The absorption line -shape may be understood within a 
phenomenological "effective liquid" approximation. The 
extra absorption producing the prominent peak in Fig. 15 
is due to the loss of a whole layer from the solid into 
the liquid (Fig. 16). For small Vo, the confined solid re- 
sponds by centre of mass fluctuations (q — ► phonons) 
shown by oscillations of N s with time (Fig. 18). Scat- 
tering from this and other sources [42-45, 48] constitute 
a background which we ignore, as a first approximation, 
for simplicity. Within our approximation, the momen- 
tum loss at the interface is modelled as regular dissipa- 
tion within a liquid strip of (fictitious) width £. The ex- 
pected momentum transfer at the interface H e = HqX the 
probability that the momentum Ho required to eject the 
layer, exists (see Fig. 19). If a local "temperature" T; oca / 
measures the degree of (de)coherence of the momentum 
transfer, then U e = (l/2)n erfc[(n - V ) / 'V2k B T local ] 
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FIG. 18: A plot of the total number of particles N s within 
the solid region (fi = 4.8) as a function of time for Vo = 1 
(top) and 6. Note oscillations in N 3 ; only the stronger pulse 
changes the number of solid layers from 22 to 21. 



FIG. 20: Average pulse velocity c pu i se (Vo) for [i — 4.8; note 
the dip in c pu i se where absorption is strongest. 
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FIG. 19: (Color online) A schematic diagram showing the mo- 
mentum transfer as assumed in our phenomenological theory. 

and £ may be extracted from Vq — H e = Vo cxp(— aw 2 £). 
Substituting for II e we obtain the extra absorption due 
to the interface, 

A 2 =4.a C ^ = -alog[l-^erfc{-^^}] (12) 

iVo \/ Iks ± local 

We use a, IIo and Ti oca \ as fitting parameters. In Fig. 15 
we show a fit to Eq. 12 of our MD data and observe that 
it reproduces all the features remarkably well. The larger 
error-bars near the peak in A 2 (Vb) reflects the difficulty 
of fitting a Gaussian to the transmitted pulse when dis- 
sipation is large. Indeed, in this region the pulse shape is 
systematically distorted away from Gaussian due to ef- 
fects beyond the scope of our simple theory. Large fluc- 
tuations in N s (Fig. 18) lead to expected [42-45] and 
detectable decrease in average pulse speed (Fig. 20). 

VI. HEAT TRANSPORT 

In this section we focus on the transport of energy 
across the system of forced solid in contact with its own 



liquid. The heat conductivity A is defined by the cele- 
brated Fourier's law j^; = — AVT where }e is the heat 
current density and VT is the temperature gradient. The 
transport of heat through small and low dimensional sys- 
tems has enormous significance in the context of design- 
ing useful nano-structures [48] . A large number of recent 
studies in lower dimensions has shown that heat conduc- 
tivity is infact divergent as a function of system size [52- 
54]. Thus it is more sensible to calculate the heat cur- 
rent or conductance of the system directly, rather than 
the heat conductivity. Therefore we focus on the heat 
current densities j& ( > ) flowing across the solid-liquid 
interface from high to low temperature and heat conduc- 
tance G = j E /AT (or resistance R = 1/G), AT ( > ) 
being the temperature difference between the two edges 
of a given region. In this study, we are particularly inter- 
ested in exploring the impact of structural changes, viz. 
the layering transitions, on heat transport. This is to re- 
member that, layering transitions are strongly dependent 
on the small system sizes and gets washed away as one 
goes to larger systems. Recently, electrical [55] and ther- 
mal transport [21] studies on confined solid strips have re- 
vealed strong signatures of structural transitions due to 
imposed external strain. Heat transport across a model 
liquid-solid interface has been studied in three dimensions 
with the interatomic potential being Lennard-Jones[25]. 
In Ref. [25] it is shown that the Kapitza resistance [26], the 
interfacial resistance, can reach appreciably large values 
when the liquid does not wet the solid. 

A. Non-equilibrium molecular dynamics 

The specific context in which we study thermal prop- 
erties of the liquid-solid interface is the same as that we 
have used for our studies of the acoustic properties in the 
last section. Again, a solid region is created within a liq- 
uid using an external chemical potential trap. We report 
results for 1200 particles interacting via the soft disk po- 
tential u(rij) = l/r^j taken within an area of 24 x 60. 
In absence of any external potential, a 2d system of soft 
disks at this density p ss 0.65 remains in the fluid phase. 
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The length scale is set by soft disk diameter d = 1, en- 
ergy scale by temperature ksT and the time scale by 
t s = y 1 md 2 /ksT. The unit of energy flux je is thus set 
by (kBT/r s d). The unit of resistance and conductance 
are r s d and (T s d) _1 respectively. The smooth interaction 
potential allows us to use standard MD simulations with 
a velocity Verlet algorithm. The time step of St = 10~ 3 r s 
in our MD ensures that the total energy is conserved (in 
equilibrium) to within 10 -4 . Periodic boundary condi- 
tions are applied in the rr-direction. We use the standard 
velocity Verlet scheme of MD with equal time update of 
time-step St, except when the particles collide with the 
'hard walled' heat reservoirs at y — and y = L y . We 
treat the collision between the particles and the reservoir 
as that between a hard disk of unit diameter colliding 
against a hard structure less wall. If the time, r c , of the 
next collision with any of the two reservoirs at either end 
is smaller than St, the usual update time step of the MD 
simulation, wc update the system with r c . During the 
collision with the walls Maxwell boundary conditions are 
imposed to simulate the velocity of an atom emerging 
out of a reservoir at temperatures T L (at y = 0) or T R 
(at y = L y )[bS\. This means that whenever a soft disk 
collides with either the left or the right wall it gets re- 
flected back into the system with a velocity chosen from 
the distribution 
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where T\y is the temperature (T& or Tr) of the wall on 
which the collision occurs. During each collision energy 
is exchanged between the system and the bath. Thus 
in our molecular dynamics simulation, the average heat 
current flowing through the system can be found easily 
by computing the net heat loss from the system to the 
two baths (say Ql and Qr respectively) during a large 
time interval r, once the system has reached steady state. 
The steady state heat current from right to left bath is 
given by (J) = lim r ^ 00 Q l /t = - ]hn T -> 00 Qr/t. In the 
steady state the heat current (the heat flux density inte- 
grated over x) is independent of y. This is a requirement 
coming from current conservation. For a homogeneous 
system ]e — (J) / L x . However if the system has inhomo- 
gencities then the flux density itself can have a spatial 
dependence and in general we can have je = 3e{x, y). In 
our simulations we have looked at Je{x, 0) and Je{x, L y ). 



B. Results 

In Fig. 21 we show that even in presence of a tem- 
perature difference across the system, the solid-liquid 
interface is formed and the liquid near the interface 
shows smectic-likc density modulation due to the pres- 
ence of nearby solid. These features are apparent from 
the structure factor calculated outside the region S (liq- 
uid), inside the region S and in the small region over 
which external potential goes to the value —fj, from zero 
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FIG. 21: (a) From left to right, structure factors in liquid, 
solid and interface regions. Interface shows clear smectic pro- 
file. Data taken at p — 13. (b) The local density profile 
along y-direction at fi — 13. (c) The isothermal compressibil- 
ity kt as a function of y at /i = 13. Compressibility shows 
strong peaks near the interfaces. Due to small size, interfacial 
enhancement of compressibility permeates right through the 
whole of solid region. In (b) and (c) lines are guides to eye. 



(Fig. 21. (a)). The local density profile also shows con- 
stant large value corresponding to the solid formed in 
region S (Fig. 21. (b)). The density of the liquid near the 
cold (ksT — 0.5) left-reservoir is higher than the density 
of liquid near the hot (fcgT = 1.5) right-reservoir. In 
Fig. 2 1(c) we plot the local compressibility Kt(v) defined 
via kt = p~ 2 (dp/d/j,)T- The compressibility of the in- 
terfaces is very large making the narrow solid region also 
unusually compressible pointing to the presence of large 
local number fluctuation. This behaviour helped in large 
shock absorption as discussed in Sec.V. Before going into 
the details of heat transport in this system, let us first 
enlist some details of the structural transitions obtained. 
With increase in the strength of the trapping potential, 
both in equilibrium and in non-equilibrium situations, we 
observe two modes of density enhancement: (a) A whole 
layer of particles enter to increase the number of lattice 
planes in y-direction. This happens, e.g., as [i is increased 
from 7 to 8. Thus in this mode the inter lattice plane sep- 
aration decreases (see Fig. 22(a)). (b) Each of the lattice 
planes grow by one atom thereby decreasing the inter- 
atomic separation within a lattice plane. This happens, 
e.g., as one increases fj, from 10 to 12 (see Fig. 22(c)). In 
the intermediate configurations one observes metastable 
dislocation pairs ( Fig. 22(d)) and peaks in the local par- 
ticle density that hops back and forth between two neigh- 
bouring positions ( Fig. 22(b)) to maintain commensura- 
bility. In the mode (a) one observes a sudden compres- 
sion in the y-direction associated with positive £d and 
7d, while in the mode (b) the system undergoes tension 




FIG. 22: (Color online) Overlapped density plot of 500 con- 
figurations in the region trapped by external potential /i: (a) 
A 23 x 23 triangular lattice solid at fi = 8. (b) Local den- 
sity peaks hop in ^-direction to incorporate > 23 particles 
in lattice planes in response to increased potential /i = 11. 
(c) A 24 x 23 triangular lattice solid at fi = 12. Notice the 
increase in particle numbers in the lattice planes, (d) Config- 
uration obtained after 150005t as a 24 x 23 steady state solid 
at /i = 16 is quenched to fi = 24. This shows a dislocation 
pair - a 23-layered region trapped in between a 24-layered 
solid. At steady state (after a time 10 St) dislocations an- 
nihilate to produce a 24 x 24 triangular lattice solid. Color 
code: blue (dark): low density and red (light): high density. 



associated with negative Sd and jd- With increase in /i, 
these two modes repeat one after another, in cycle. This 
allows the system to release the extra stress developed 
in one direction due to particle inclusion in the previous 
cycle by inclusion of particle in the other direction in the 
next one. Certainly, at large enough fi the solid region 
goes towards very high packing fraction, thereafter stop- 
ping the process of particle accumulation. To summarize 
the major structural changes obtained in soft disk solid, 
we find, strained triangular solids with 23 x 23, 24 x 23, 
24 x 24 unit cells at n = 8, 12 and 24 respectively (See 
Fig. 22). The large elastic and core energy cost inhibits 
the formation of an equilibrium nano solid with dislo- 
cations - even if dislocations form, the solid eventually 
gets rid of them by either incorporating particles from or 
ejecting them to the liquid part. 

Before any measurement is done, the system is allowed 
to reach the steady state where the current density in- 
tegrated over the whole x-range is the same at all y. 
If LTE is maintained, v(y) is expected to obey Gaus- 
sian distribution locally. That gives definitions of local 
temperature from all the even moments of v(y). Thus 
k B T{y) = (1/2 mv 2 {y)) and k B T{y) = m^/(v 4 (y)}/8 etc. 
To check for the local thermal equilibrium (LTE) from 
our simulation, we find (v 2 (y)) and (v 4 (y)) as a func- 
tion of distance y from cold to hot reservoir and compare 
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FIG. 23: (Color online) (a) Plot of temperature profile 
(1/2 mv 2 (y)) and m % /(« 4 (y))/8 at fi = 8. (b) The tem- 
perature profile (1/2 mv 2 (y)) as a function of y, the system 
coordinate perpendicular to the reservoirs, for fi = 8, 16, 21. 

the above mentioned definitions of fcsT(y)(Fig. 23(a)). 
From Fig. 23(a) it is evident that the temperature pro- 
file is almost linear in the single phase regions, like the 
liquid and the solid, with sharp increase near the inter- 
faces and the LTE is approximately valid in all regions. 
In Fig. 23(b) we plot temperature profiles ksT(y) as ob- 
tained from (l/2mv 2 (y)} at well separated trapping po- 
tentials /i = 8, 16, 21. With increased trapping strength, 
the temperature difference between the edges of the solid 
region decreases indicating an enhancement of heat con- 
ductance within the solid. The temparature jumps at 
the interfaces also increase with increasing trapping po- 
tential. Such a jump in the temperature is known as the 
Kapitza or contact resistance (Rk) [26]. This is defined 
as, 

AT 

R K = — (14) 

JE 

where AT is the difference in temperature across the in- 
terface. It is evident that the interfaces are the regions 
of the highest resistance in the system. This large resi- 
tance can be traced back to large density mismatch at the 
contact of two phases. The conductance of the high tem- 
perature liquid near the right reservoir with k B Tn = 1.5 
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FIG. 24: Plot of the heat flux, je, as a function of the trap 
depth, p. Note that the overall flux decreases as a function 
of p. is expressed in units of ksT '/r s d. 

is lower than the low temperature liquid in the other 
side in contact with the reservoir with UbTl = 0.5. The 
temperatures are expressed in units of ksT. The heat 
current, as expected, flows from the right reservoir to 
the left reservoir. 

In Fig. 24 we have plotted the heat flux through the 
system as a function of p. As p increases, the atoms 
from the surrounding liquid get attracted into the po- 
tential well and the density of the solid progressively be- 
comes higher at the cost of the liquid. Because of Kir- 
choff 's law, the heat conductance of a composite system 
of liquid-solid-liquid connected in series through interfa- 
cial regions is dictated by the low conductance regions. 
The decrease in liquid density decreses the pressure in the 
liquid regions, thereby reducing the heat conductance in 
them [56]. The conductance in liquid is always lower than 
solid. These result in an overall decrease in the heat flux 
and consequently the overall conductance. Moreover, at 
larger p there occurs larger density mismatch at the inter- 
faces leading to larger Kapitza resistance (Fig. 26). The 
change in density is sharper near the two layering tran- 
sitions - thus heat flux shows sharper drops near the 
transitions at p = 7 — 8 and p = 11 — 12. However, 
close to the layering transition at p <~ 8 there is a local 
peak in the value of the heat flux suggesting that a sig- 
nificant amount of kinetic energy is exchanged between 
the liquid and solid through the interface at the layering 
transition. This excess conduction is due to an enhanced 
number fluctuation in the direction of heat flow in this 
mode of layering transition. 

In Fig. 25 we show the heat conductance in solid region 
G s as a function of strength of the trapping potential 
p. The inset in Fig. 25 shows the change in the averaged 
density of the solid region p s d 2 . The p s d 2 — p plot show 
clear staircase-like sharp increases near p, = 8, 12. As p 
is increased from 7 to 8, a layering transition in a direc- 
tion perpendicular to the interfaces occurs; whereas as 
p is increased from 10 to 12 each of the lattice planes 
lieing parallel to the solid-liquid interfaces grow by one 
atom (see also Fig. 22(a) & (c)). These two modes of den- 
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FIG. 25: Plot of the thermal conductance of the solid region, 
G s as a function of p. G s is plotted in units of (r s ci) _1 . The 
inset shows change in solid density p s d 2 as a function of p. 
Jump increase in p s d associated with layering shows up in 
sudden large increase in conductance near p = 7.5. Another 
sharp increase in p s d 2 near p = 12 is due to growth of lat- 
tice planes by one atom each; this happens in the orthogonal 
direction to heat conduction and does not affect G s - 



sity fluctuations leave their signatures by enhancing heat 
conductance G s . Notice that, the layering transition at 
p = 8 increases stress in the solid region in the direction 
of heat conduction, thereby showing a step-like increase. 
The other mode of density fluctuation is in the normal 
direction to heat transport and thereby affects the heat 
conductance only near the transition due to the associ- 
ated enhancement of overall fluctuations (a sharp peak 
at p = 12). 

We now find out the Kapitza resistance across the solid 
liquid interface as a function of the strength of the exter- 
nal potential p. With increase in p, the system shows 
a jump in the density of the solid region correspond- 
ing to the addition of an entire layer of atoms (see in- 
set of Fig. 25). From the profile shown in Fig. 23(b), 
the Kapitza resistance is easily obtained by dividing the 
temperature jump by the energy flux. Note that a slight 
dependence of Rk on T is visible in Fig. 23(b) with a 
larger temperature jump on the "warm" side. The re- 
sults shown in Fig. 26 correspond to the average values 
of Rk over the "warm" and "cold" sides. The plot of 
Rk as a function of p shows a distinct jump as a layer is 
included (for a value of p close to 8) in the solid region 
(Fig. 26). The jump in Rk is also accompanied by a 
local dip at the transition, corresponding to larger num- 
ber fluctuations at the interface. The combined effect of 
the enhanced Kapitza resistance as well as enhanced con- 
ductance of the solid can be summarized by defining the 
Kapitza length in units of width of the solid region L s as 
Ik = RkG s . This is a measure of excess width of a solid 
which is equivalent in giving rise to a resistance equal to 
the Kapitza resistance. This is reminiscent of the "effec- 
tive liquid" concept which has been used in explaining 
some of the features of acoustic shock absorption in the 
last section. The Kapitza length also shows a peak near 
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region within this time scale when the solid is decorated 
by the dislocation pair. This gives a heat conductance 
G s = 2.29 (r s d) _1 . After a further wait for 10 5 St the dis- 
locations get annihilated. At this stage the whole trapped 
region is transformed into a 24-layered solid. Then the 
heat conductance comes out to be G s = 3.53 (r s d) _1 . 
Thus after complete annihilation of the dislocation pair 
the conductance of the solid rises by about 54%! This 
study already indicates towards the fact that dislocations 
behave like large resistances towards heat transport. One 
can use the presence of stable dislocation-pairs as have 
been obtained in confined narrow strips [20, 21] for a more 
careful study of impact of dislocations on heat transport. 



FIG. 26: Plot of the Kapitza resistance, Rk, expressed in 
units of r s d as a function of fx, shows a jump at the layering 
trans 
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C. Approximate theory 

We now provide an approximate theoretical approach 
to calculate heat conductance within the solid region. 
We use a free volume type calculation to obtain an ap- 
proximate estimate for heat conductance starting from 
an exact expression for a-th component of the heat flux 
density[21] 

j a (r) = j*(r)+#(r) = £j(r- ri )fcv? 

% 

(15) 



FIG. 27: Plot of Kapitza length Ik in units of L s , as a function 
of fj,. This shows a jump increase at the layering transition. 



the layering transition at /i = 12 (Fig. 27). 

The layering transitions in solid occur via dislocation 
formation which ultimately annihilates by incorporat- 
ing more particles from the liquid region. The kine- 
matics of dislocation formation and annihilation is as- 
sisted by diffusion and dislocation climb which are very 
slow processes [20] in a solid compared to particle colli- 
sion and kinetic energy transfer times. Thus it is pos- 
sible for a system with metastable dislocation pairs to 
reach a thermal quasi-steady state. Fig. 22. (d) shows 
overlapped configurations of the solid region contain- 
ing a dislocation-antidislocation pair, as the system is 
quenched from \i = 16 to fi = 24. The overlapped config- 
urations are separated by time 100 St and collected after 
a time of 15000 St after the quench begins. Since dislo- 
cations annihilate though a conserved diffusive dynamics 
which takes a long time (10 5 St) compared to the particle 
collision and kinetic energy transfer times, the system 
in presence of dislocation pairs is in a effective steady 
state. It also maintains LTE that we have checked by 
computing (v 4 (y)} and (v 2 (y)} locally. Thus, in a similar 
manner as stated before, we find the conductance in solid 



obtained from continuity of local energy density. Here 
6(x) is the Heaviside step function and S(. . . ) is a Dirac 
delta function, hi = 77w t 2 /2 + </>(r.;)+ J2i>j u ( r ij), 4>{ v i) is 
an onsite potential and u{rij) is the inter-particle inter- 
action. The first term in Eq.15 j& ( r ) denotes convection 
while the second term j%(r) denotes conduction. The 
above formula for conduction has an simple interpreta- 
tion. The sum is over only those i for which xf > x a . 
Thus this formula gives the net rate at which work is 
done by particles on the left of x a on the particles on the 
right which is therefore the rate at which energy flows 
from left to right. The a-th component of the integrated 
heat current density over the solid region [21], 



du{ rij ) x%x% g | ^ 



(16) 



In solid region most of the transport is carried out 
by conduction and one may ignore the convection part. 
In this study we focus on the average heat current den- 
sity along y-direction, je = (I y ) / L X L S . Ignoring convec- 
tion inside solid, approximating the system as a system 
of hard disks with some effective hard disk diameter a 
and assuming LTE [21], the heat conductance in units of 
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(r s d) 1 can be expressed as, 



G s 



JE 

AT 



,Ps Vc 



(17) 



where p s = 4ry s /7r, y c is the average separation between 
the colliding particles in y-direction, t c is the mean col- 
lision time and p s is the average density of the solid. For 
the details of this derivation refer to Ref. [21]. The extra 
factor of (d/a) 2 is due to the mapping of the soft disks 
of diameter d to effective hard disks of diameter a. 

Now this conductance G s can be calculated if one can 
obtain some estimate for y c and r c . We estimate y 2 and r c 
from the fixed neighbor free volume theory (FNFVT) as 
in Ref. [21]. For different values of (rj s ,L s ) [ i] s obtained 
by extremizing the free energy of the full system at a 
given potential well depth /j,, as discussed in Sec. Ill ], one 
can obtain ao, e xx , e yy that gives the basic geometrical 
inputs of b, and h (Fig. 9). We assume the test particle Pq 
moves in the cage formed by its neighbors and obtain the 
average values [y 2 ]f v and [r c ]/„ for the moving particle 
from FNFVT. We assume that the position of the centre 
of the moving disk Po(x, y)i at the time of collision with 
the other disks, is uniformly distributed on the boundary 
B of the free volume. Then [y 2 ]f v can be easily calculated 
using the expression [21] 



Lb 



(18) 



where Bi is the part of the boundary B of the free volume 
when the middle disk is in contact with the i th fixed 
disk, ds is the infinitesimal length element on B while Lb 
is the total length of B. An exact calculation of [r c ]/„ 

is nontrivial. However we expect [t c ]/„ = c v 1 ^ 2 /T 1 / 2 
where Vf v is the "free volume" [see Fig. 9] and c is a 
constant factor of 0(\) which may be used as a fitting 
parameter. Thus the heat conductance in solid region 
may be expressed as, 



[G s ]f v = 



3 Ps a 2 T^ 2 [y 2 



fv 



LA 2 



c v 



1/2 
fv 



(19) 



There are well defined schcmes[57] to approximate a soft 
disk system as a system of effective hard disks. How- 
ever, we choose a much simpler path of finding the heat 
conductance of a trapped hard disk system instead. In 
Fig. 28 we plot [G s ]fv as a function of \i with L s = 19.73cr 
the same width of the trapped hard disks used in Sec. II. 
For hard disks we obtain heat conductance in units of 
(rer) -1 and use an average temperature UbT = 1 and 
set c = 1. This gives the estimate of heat conductance 
along y-direction in a hard disk system that showed a 
layering transition from 21 to 22 layers near \i — 8 in 
Fig. 7. The plot clearly shows the associated increase in 
heat conductance in the solid region. This behavior is 
also in qualitative agreement with Fig. 25. 

With the help of these results we may conclude that 
the layering transition has a profound effect on the ther- 
mal properties of the trapped solid licing in contact with 
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FIG. 28: Free volume estimate of the y- component of heat 
conductance [G s ]j v , in units of in the solid region 

of a hard disk system composed of liquid and trapped solid 
region. Heat conduactance shows jump increase as layering 
transition occurs with increase in trapping potential p. 



its liquid. An important consequence of this study is the 
possibility that the thermal resistance of interfaces may 
be altered using external potential which cause layering 
transitions in a trapped nano solid. Moreover, the heat 
conduactance of the solid may be drastically reduced by 
tuning the trapping potential. We believe that these phe- 
nomena have the potential for useful applications for e.g. 
as tunable thermal switches or in other nano engineered 
devices. 



VII. CONCLUSION 

In this paper we have shown that liquid-solid interfaces 
which are constrained by strong chemical potential gra- 
dients remain flat and at the same time undergo fluctua- 
tions which increases or decreses the number of solid lay- 
ers by one. The nature of these fluctuations are strongly 
influenced by the size of the solid. It is expected that 
for macroscopically large solids these fluctuations would 
reduce to the random nuclcation of steps on the solid 
surface and the sort of coherence observed here would be 
absent. Confining a thin 'long' strip of solid by smooth 
walls in quasi one dimension leads, strictly speaking, to 
a destruction of solid- like order [20, 35] and an enhance- 
ment of smectic-like ordering of individual layers parallel 
to the confining walls. It is this reduction of interlayer 
coupling which is ultimately responsible for the spallation 
of single solid layers. 

Before we end, we would like to examine carefully the 
relevance of our results to practical situations and ex- 
perimental systems constructed in a laboratory. In re- 
cent years our ability to manipulate matter at an atomic 
or molecular level has increased tremendously [58]. It is 
possible now to localize atoms using carefully designed 
atomic traps[30, 59] and observe their properties. It is 
also possible to set up experiments where atoms may 
be picked one by one and arranged in any specified 
pattern [58]. Parallel to this development, one can now 
synthesize functionalized colloidal[27, 60] particles with a 
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variety of shapes and sizes which mimic the properties of 
atomic matter at length and time scales which are easy 
to handle even in relatively inexpensive experimental set- 
ups. Colloidal particles can also be manipulated using 
laser tweezers and traps [19, 29, 30, 32], confining within 
narrow channels and slits[19, 61], adsorbing on substrates 
or air- water interfaces [62] and assembling layer by layer 
using carefully designed substrate templates[28, 63]. 

The structural aspects of our results viz. the layering 
transition[18, 19] and all associated phenomena should 
be observable in both atomic systems like rare gases in 
atom traps and in colloids on templates or in laser fields. 
Indeed, interfacial fluctuations similar to the sort dis- 
cussed in this paper have been observed during early 
nano-indentation experiments [64]. Layering transitions 
have also been observed in shaken hard disks more than 
two decades ago [65]. In confined molecular systems such 
layering transitions are of great relevance to the study of 
nano-tribology[18] . 

The dynamical aspects of our results, on the other 
hand, will be difficult to observe in colloidal systems be- 
cause of viscous damping by the solvent. Nevertheless, 
it may still be possible to observe some of these effects 



if this damping is small[66]. Our results for momentum 
and heat conduction in trapped solids therefore pertain 
mainly to atomic or molecular systems where such damp- 
ing is absent [18]. As such, we do not see any difficulty for 
generalizing the main conclusions of our study to dimen- 
sions higher than two. In three dimensions, confinement 
to a slit should have similar effects on a three dimensional 
solid viz. reducing inter Iyer coupling so that two dimen- 
sional layers may become partially independent. Weak 
shock waves may then be able to spallate these layers 
into a surrounding liquid or gaseous phase. Work along 
these lines is in progress. 
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